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Abstract 

In this paper we solve numerically the two dimensional elliptic sine-Gordon equation with appro- 
priate boundary conditions. These boundary conditions are chosen to correspond to the Josephson 
interaction between two adjacent pancakes belonging to the same flux-line in a highly anisotropic 
superconductor. An extrapolation is obtained between the regimes of low and high separation of 
the pancakes. The resulting formula is a better candidate for use in numerical simulations than 
previously derived formulas. 
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I. INTRODUCTION 



High-temperature superconductors belong to the class of superconducting materials 
known as type II that allow for partial magnetic flux penetration whenever the external 
field satisfies H cl < H < H ^ 1 ^ 2 ^ . The flux penetrates the sample in the form of flux-lines 
(FL's), each containing a quantum unit 0o — hc/2e of flux. At low temperature the FL's 
form an ordered hexagonal lattice (Abrikosov lattice) due to their their mutual repulsion. At 
high temperature and/or magnetic field this lattice melts due to thermal fluctuation o 4 i 5 i 6 i 7 i 8 . 

High-temperature superconductors are anisotropic material which are made from stacks 
of superconducting layers associated with Cu0 2 planes. The layers are weekly coupled to 
each other. The parameter measuring the anisotropy is 7, defined as 7 2 = m z /m±, where m z 
and m± denote the effective masses of electrons moving along the c axis (perpendicular to the 
superconducting planes) and the ab plane, respectively. While for the material YBa 2 Cu30 7 _<5 
known as YBCO the anisotropy is somewhere between 5-7, for the material Bi 2 Sr 2 CaCu 2 0g +( 5 
known as BSCCO, the anisotropy is estimated to be between 10 to a 100 times larger. 

For BSCCO and highly anisotropic materials similar to it, each FL is represented more 
faithfully by a collection of objects referred to as "pancakes". Pancakes are centered at 
the superconducting planes. Each pancake interacts with every other pancake, both in the 
same plane and in different planes. The interaction can be shown to consist of two parts. 
The first part is called the electromagnetic interaction (or simply magnetic) and it exists 
even in the case that the layers of the materials are completely decoupled, so no current 
can flow along the c-axis of the sample. A pancake vortex located in one plane gives rise to 
screening currents in the same plane as well as in all other planes. A second pancake vortex, 
located elsewhere, interacts with the screening currents induced by the first pancake^. This 
interaction has been calculated by Clem and others^. Two pancakes in the same plane 
interact with a repulsive interaction while pancakes in different planes attract one another. 

The second part of the interaction among pancake vortices is the so-called Josephson 
interactio n 2 ! 9 ! 12 . It results from the fact that there is a Josephson current flowing between 
two superconductors separated by an insulator and this current is proportional to the sine 
of the phase difference of the superconducting wave functions. The superconductors in the 
present case are the different Cu0 2 planes. When two pancakes belonging to the same 
stack and residing in adjacent planes move away from each other, the phase difference that 
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originates causes a Josephson current to begin flowing between the planes. This results in 
an attractive interaction between pancakes that for distances small compared to r g = jd is 
approximately quadratic^ in the distance. When the two adjacent pancakes are separated 
by a distance larger than r g , a "Josephson string" is formed, whose energy is proportional 
to its length 12 . 

When doing simulations of flux lattices it is important to include correctly the strength of 
the interaction among pancakes. It turns out that for anisotropies smaller than about 150, 
the Josephson interaction is dominant over the magnetic interaction and the later can be 
included only in the form of an effective in-plane interactio n 2 ' 13 ' 15 . For higher anisotropies 
the magnetic interaction must be included among all pairs of pancakes , but the Josephson 
interaction still matters for any finite anisotropyJ^. The form of the Josephson interaction 
used in simulation was first introduced by Ryu et alM based on the Lawrence- Doniach 
model 16 . These authors used a certain approximation which can be somewhat improved. It 
is the aim of this paper to review the approximations made for the Josephson interaction 
and to suggest a better approximation to be used in numerical simulations. In order for 
future Monte Carlo simulations to yield a better agreement with experimental results it is 
important to choose the form of the interaction among pancakes to be as close as possible to 
the exact interaction in real materials. In the simulations the interaction among pancakes 
belonging to different FL's is dominantly electromagnetic for the range of magnetic fields 
used in experiments in the vicinity of the melting transition. For the same FL, the Joseph- 
son interaction is present predominantly for nearest neighbor pancakes which are displaced 
from an alignment along the z-axis. For pancakes separated by more than one plane the 
electromagnetic interaction is present. If a kink is present in more than one place along 
the same FL the total interaction is given to a good approximation as a sum of the pair 
interactions, provided the FL does not deviate too much from a straight line which is usually 
the case on the solid side of the melting line. 
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II. THE MODEL 



Our starting point is the Lawrence-Doniacb2J£ Gibbs free-energy functional, 
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where ip n represents the superconducting order parameter in the n th Cu02 layer, is 
the vector potential in the plane, and a z its component perpendicular to the planes, b is 
the local magnetic field and H is the externally applied field, a and /3 are the Ginzburg- 
Landau parameters, m and M are the effective masses along the x — y and z directions 
respectively, d is the separation between the superconducting planes. <fio is the flux quantum. 
The integration along the z-direction has been replaced by a discrete summation over the 
superconducting layers. In the London approximation the absolute value of ipn is treated 
as a constant and its phase (p n is varying. The Gibbs functional becomes (dropping some 
constants): 
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where A is the penetration depth and and 7 = \jM/m is the anisotropy. Minimization with 
respect to a^ 2 ^ gives 
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where A is the 3-dimensional Laplacian. Minimization with respect to a z gives 
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is the gauge invariant phase difference between the layers n and n + 1, and 
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is the Josephson-coupling current density between layers. Minimization with respect to n 
gives 

A®</> n + ^V (2) • a (2) = M*»,«-i) - sin($ n+1>n )] . (2.8) 

We have used the Coulomb gauge V • a = 0. Eqs. ()2.4j) . (j2.5J) and ()2.8|) are to be solved 
with the appropriate boundary conditions and the solution must be substituted back into 
the expression for the Gibbs free-energy. 

An isolated pancake residing in plane n is a singular solution of the equation for the phase 
of the wave function which -n the limit of R — > R n satisfies 



V( 2 WR) = - gx(R R " 



. (2.9) 
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where R n denotes the center of the pancake and z is a unit vector in the z-direction. This 
solution becomes exact in the whole plane in the infinite anisotropy limit. As one fully 
encircles the pancake the phase n changes by 2tt, and is singular at the center of the 
pancake. In the case of infinite anisotropy, the full solution of Eqs. ()2.4j) . (|2.5jl and ()2.8|) 
can be founo-^ and from it one can deduce the interaction energy for two pancakes in the 
same plane or in different planes. But for more complicated configurations, or when the 
anisotropy is finite, the solution can only be found approximately. 

Consider Eq. ()2.8|) for the (n + l)-layer and for the n'th layer respectively. Subtracting 
the second case from the first we obtain 
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adding and subtracting the term {2nd/ ipo) Aa z it becomes 
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using the coulomb gauge we see that 
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We are now in a position to simplify Eq. (|2.11j) by using Eq. (|2.6|) and Eq. (|2.5j) together 
with Eq. to finally fincU^ 

A (2) $ n+ i, n = [2sin($„ + i >n ) - sin($ n+2 ,n+i) - sin($ n , n _i)] 



which amazingly involves only the gauge invariant phase difference among the layers. So far 
this equation is nearly exact. The problem is that it represents an infinite set of coupled 
equations. 

Consider now a FL which consists of a stack of pancakes. Assume that there is a kink in 
the FL in the sense that the n'th and the (n + l)'th pancakes are not on top of each other 
but are shifted a distance w along the x direction. The pancakes with index > n + 1 are 
all aligned along the ^-direction as well as the pancakes with index < n. See Figure (JJ). 
In this situation one might think that except for $ n+ i )n all the phase differences vanishes. 
However this is not entirely true since because of the Josephson coupling there is induced 
phase differences away from the location of the kink. For example in the situation that 
w is very large a so called Josephson vortex is formed. Although its core only extends a 
distance d in the ^-direction and a distance 70? in the ^-direction, the associated magnetic 
field extends a distance A in the ^-direction and 7A in the y-direction. 

Nevertheless in order to truncate the infinite set of equations we will assume that we 
can neglect all the gauge invariant phases as compared to $ n+ i jn . This turns out to be 
quite satisfactory for distances much smaller than 7A. Later we will add an overall factor 
to partially compensate for the approximation made. If we make this approximation, and 
further neglect the last term on the right hand side of Eq. ()2.13|) since A 2 3> d 2 , we obtain^ 



which is the well-known elliptic sine-Gordon equation in two dimensions, where we put 
$ = $ n+ i in and A = 7c/. The energy associated with a solution to this equation due to the 
Josephson currents is found by evaluating the expression 
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FIG. 1: A kink in a configuration of pancakes aligned along the z-direction. 



The boundary conditions for the sine-Gordon equation appropriate to the configuration 
under discussion are fixed on the geometry depicted in Fig. (J2J). The outer boundary is a 
large circle of radius R m which tends to infinity. On this boundary we set $ = 0. The inner 
boundary consists of two tiny circles the centers of which are a distance w apart, connected 
by a thin rectangle. On the inner boundary we take the solution to satisfy 



where R = (x, y) and atan2 is the 4-quadrant inverse tangent whose values lie in the interval 
(— 7r, 7r) as depicted in Fig. (j3J). It is to be distinguished from arctan(?//x) whose value ie in 
the interval (— 7r/2, 7r/2). Thus the values of $ is given approximately by the values depicted 
in Fig. These values become exact as the radius of the inner circles and the width of 
the rectangle tend of zero. 

We proceed by first discussing the solution under certain limiting conditions. Under these 
conditions we can come up with analytic expressions. One such limit is the case w < A. 
Under this condition it is a good approximation to linearize the sine-Gordon equation and 
solve instead the equation 



$(R) = atan2(y, x — w/2) — atan2(y, x + w/2), 



(2.16) 



(2.17) 



which has a solution 




(2.18) 
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FIG. 2: Boundaries and boundary conditions for the solution of the two-dimensional elliptic sine- 
Gordon equation. 




FIG. 3: definition of atan2 



where (R, 8) are the polar coordinates of the vector R and K\ is the modified Bessel function 
of order 1. With the given boundary condition this solution is only valid for R 3> w. We 
see that for w <C R <C A the solution becomes 



$(R 

whereas for i?>A we obtain 
$(R) - 



~ — sin(0), 



R 



w <C R < A, 



(2.19) 



To obtain the energy we note that since $ is small we can expand the cosine in Eq. (|2.15J1 



to second order in $ and substituting from Eq. (|2.19|) cutting the integration at w for lower 
limit and A for an upper limit we obtain for the potential energy as a function of w 

The reader may wonder how the normalization of the solution in Eq. (J2.18j) has been de- 
termined. This is due to the boundary conditions on the inner boundary. In the region 
u» < E < A the term on the right hand side of Eq. (|2.17j) can be dropped and <3> satisfies 
Laplace's equation. A solution that satisfies the boundary conditions can be written in the 
first quadrant in terms of the ordinary arctan function as follows: 

$(R) = arctan ( V —— J- arctan ( V —— } , (2.22) 

\x — w/2J \x + w/2J 



which for x ^> w gives 



*(R)*^^sin W , (2.23) 



in agreement with Eq. (|2.19|) . The solution given in Eq. (|2.18|) matches correctly to the 
solution of the Laplace's equation satisfying the correct boundary conditions and hence it 
is properly normalized. 

The other asymptotic limit that can be done analytically is the limit w — > oo. In that 
limit the solution should be independent of x, and thus $ should satisfy the one-dimensional 
sine- Gordon equation 

^a*(y) = ^Bm($(y)), (2.24) 
with the solution in the upper half plane satisfying 

$(0+) = vr, $(oo) = 0, (2.25) 

and on the lower plane 

$(0") = -7T, $(oo) = 0. (2.26) 
A well-known kink solution to the sine-Gorgon equation is given by 

$(y) = sign(y) x 4arctan(exp(-|y|v / 2/A)), (2.27) 

which implies 

sin($(y)) = 2 sinh^v^/A)/ cosh^yv^/A). (2.28) 



The energy of this solution per unit length in the x direction is given by 

dy (l-cos(4arctan(exp(-|y|v / 2/A)))) = ^|^- (2-29) 

If w is large but finite, we can take the total length of the "string" as w, and to a good 
approximation 

"<»> = ^(i)- ra>>A - (2 - 30) 

which is linear in the separation between the two singularities (pancakes). In this regime 
the non-linearity of the sine-Gordon equation is crucial. In the case of w — > oo, which can 
be referred to as a single Josephson vortex, it is possible to find an approximate solution 
directly to the infinite set of equations given by Eq. ([2.130 . The solution found by Clem, 
Coffey and Hao 11 ' 12 reads 
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where we put 



A c = 7 A, (2.33) 



The value m = corresponds to the position of the kink (Josephson vortex). Note that 
in the ^-direction the phase difference decreases rapidly but not abruptly. Also there is an 
additional length scale A c ^> A not present in the approximation made above. If we plot the 
rhs of Eq. ([2.31)1 for m = on the same graph as sin($) from Eq. ([2.28J1 for say A = 1 and 
A c = 100 we see, as depicted in Fig. that both rise linearly in y up to a maximum for 
y ~ A and then both decay exponentially, on scales A c and A respectively. When evaluating 
the energy with the solution given by Eq. ([2.3 1[) one finds instead of the result represented 
by the rhs of Eq. ([2.29)1 a value of 

In (£)+ 1.12)^, (2.34) 

per unit length of the vortex in the x-direction. Using a different approximation Koshelevi^ 
claims that the constant in the last equation should be 1.55 instead of 1.12. In any case the 

10 



2 4 6 8 10 

Y 



FIG. 4: Comparison of the solution given by Eq. (|2,31|) (longer tail) to the solution given by 
Eq. (|2,28|) plotted against y in units for which A = 1 and A c = 100. 

appearance of \n(X/d) is the result of the fact that the magnetic field decays over a distance 
of A and not d in the ^-direction as was implied by the neglecting the phase differences away 
from the $ n +i,n- So we can correct this coefficient by introducing the logarithm by hand in 
front of the solution of the sine-Gordon equation as will be discussed further below. 

III. A NUMERICAL SOLUTION 

The main problem is how to interpolate the potential energy of two pancakes as their 
separation changes from to <C A to «) > A. In this case an analytical solution is not 
available and one has to seek a numerical solution. This we achieved using the finite-element 
method implemented by the PDE tool of the program MATLAB. We solve the sine-Gordon 
equation 12.141 m the plane with the boundary condition displayed in Fig. (J2J) and on the 
inner boundary given more accurately by Eq. (|2.16J) . We have chosen A = 1. The outer 
radius has been taken to be R ma x = 15 for 0.2 < w < 10, and R ma x = 5 for w < 0.2. The 
radius of the inner circles has been taken to be R m i n — 0.04 for 1.2 < w < 10, R m in = 0.02 
for 0.6 < w < 1.2 and R m in = w/30 for w < 0.6. The width of the rectangle has always 
been taken to be R min /2. 

In Figure El we show example of the triangular grid used for the case w — 10 and the 
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FIG. 5: Triangular grid for the case uu = 10. Distances are in units of A. 

corresponding solution &(x,y) is depicted in Figure El This figure depicts contour lines of 
equal $ for some specified values as indicated in the figure caption. It can be verified that 
the contour lines for the range of x in the middle between the singularities (—3 < x < 3), are 
almost identical to those obtained from the sine-Gordon kink solution given by Eq. (|2.27J1 
which is valid for w — > oo. For this solution the contour lines are always parallel to the 

X~clXlS. 

The energy versus pancake separation w/A is plotted in Figure A good fit to the 
energy plot is as follows: 

Mv)-<)M")< 

= (f) (2.828 (f) -1.414), 2A< m , (3.1) 

Notice that the factor 2.828 = A/y/2 agrees with the coefficient of the analytical solution 
given by Eq. ()2.3()j) . The prefactor of the quadratic term is smaller than that given by 
Eq. (j2.21|) since we made a fit up to a value of w larger than A in order to obtain a simpler 
extrapolation formula. The constant inside the logarithm captures the quadratic term on 
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FIG. 6: Contour plot of the solution of the sine-Gordon equation for w = 10. Contour values from 
bottom to top: O(boundary), -0.01, -0.05, -0.1, -0.5, -1, -1.5, -2, -2.5, -3, 3, 2.5, 2, 1.5, 1, 0.5, 0.1, 
0.05, 0.01, (boundary). 

the right of Eq. (|2.21j) which is 0((w/A) 2 ) and whose coefficient is not determined by the 
simple argument given. ( For very small values of w a good fit is given by {wV/eod) = 
1.57 (w/A) 2 ln(2.6A/w) ). In order to compensate for the approximation made by the use 
of the sine-Gordon equation, we will now rescale the energy by a factor that will reproduce 
correctly the Josephson vortex energy as given by Eq. (|2.34|) which takes the exponential 
decay along the ^-direction into account. Thus we suggest the following expression for the 
Josephson energy among two adjacent pancakes belonging to the same FL and displaced 
horizontally by an amount w: 

V new (w) = e d (1 + ln(A/d)) 0.25 (w/A) 2 In (9A/w) , w < 2A 

= e d (l+ln(A/d)) {{w/A) -0.5), 2A < w. (3.2) 

This expression should be compared with the extrapolation proposed by Ryu et a/.—, that 
reads, after correcting for a misprint by a factor of 2/ir in that reference, and shifting by a 
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FIG. 7: The energy as obtained from the numerical solution of the sine-Gordon equation (inverted 
triangles) and a fit (solid line) as given by Eq. 



constant 



V ryu (w) = e d (1 + ln(A/d)) 0.25 (w/A) 2 
= e d (l + ln(A/d)) (KA)-l) 



w < 2A 
2A < io. 



(3-3) 



These authors did not take into account the logarithmic dependence ln(A/u>) of the 
potential as given by Eq. (|2.21j) and kept only the quadratic dependence on w. They 
adjusted the coefficient to obtain a match of the function and its first derivative at the 
matching point. Note that our formula Eq. ()3.2|) also has the feature of a continuous first 
derivative (related to the force between pancakes) at the matching point. For BSCCO the 
value of ln(X/d) ~ 4.7. Thus the additive factor of 1 in the coefficient is rather small 
compared to the logarithm and this value was used by Ryu et al. instead of the more 
precise (but still approximate) value 1.12 derived in Ref. [12J. Koshelevi^, using a different 
approximation, claims that this value should actually be 1.55 and one might consider using 
this value instead of the additive 1 in the suggested new formula. 

In Figure |H1 we show a comparison of the two formulas given by Eq. ()3.2|) and Eq. (j3.3J) , 
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FIG. 8: Comparison of the two formulas given by Eq. (|3.2|) and Eq. (j.3,3j) for V new (upper curve) 
and V ryu (lower curve) with the prefactor omitted. A was taken to be equal to 1. 

where the coefficients €q(1(1 + ln(X/d)) have been omitted. The new interpolation should 
provide a better fit to the true potential in real systems. In the next section we discuss the 
results of simulations we have carried out with the new formula versus the old one. 

IV. NUMERICAL SIMULATIONS USING THE NEW INTERPOLATION 

We have carried out numerical simulations using the new formula for the Josephson 
interaction between pancakes. We also included the electromagnetic interaction along the 
lines discussed in our previous publication 15 . This simulation has been done with 36 planes 
and 36 pancakes per plane. The temperature dependence of the penetration depth used for 
the results displayed below is A 2 (0)/A 2 (T) = 1 - T/T c with A(0) = 1700 A, T c = 90 K and 
d = 15 A . The simulation method used is a multilevel Monte Carlo, as discussed in our 
previous publication a 14 ' 15 , with the only difference being the formula used for the Josephson 
interaction. 

In Figure EH we show the results of a simulation for BSCCO using the parameter 7 = 250 
for the anisotropy and B = 100 Gauss for the applied field (curves labeled as 'new'). This 
is to be compared with the results using the old formula for values of 7 = 166 and 7 = 250 
(curves labeled as 'ryu'). It should be emphasized that the figures refer to as 'ryu' are not 
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results of Ryu et. at^ but rather our simulations using their formula (corrected by a factor 
2/71") for the Josephson interaction. The figures show the structure factor and the energy 
as a function of temperature. The absolute value of the energy is not meaningful since a 
constant has been subtracted, only energy differences are meaningful. 

We see that roughly speaking the location of the melting transition for 7 = 250 using the 
new formula occurs at the same temperature as for 7 = 166 with the old formula, although 
some details like the magnitude of the structure factor and the magnitude of the energy 
jump are not identical so it not just a trivial shift. The melting transition with the new 
formula appears a little sharper. But roughly speaking there is a scaling of the anisotropy 
by a factor ~ 1.5 compared with using the old formula. Thus our previous estimates for the 
anisotropy of experimental samples should be roughly increased by this factor. In Ref. [15 
we estimated the anisotropy of experimental samples used in the experiments of Majer 
et alM to be in the range 250-450 depending on the actual temperature dependence of the 
penetration depth. This value should thus be increased to the range 375-675 with an average 
of 525. This is in good agreement with recent measurements that places the anisotropy in 
experimental samples to be about 550^. All the other qualitative results of our previous 
simulations^ should remain roughly the same apart from the fact the curves should be 
relabeled to correspond to a higher anisotropy. 



V. CONCLUSIONS 



In this paper we derived a new approximate formula for the Josephson interaction between 
neighboring pancakes belonging to the same flux-line. The formula extrapolates better 
between the limiting analytic solutions of the sine-Gordon equation than the previously 
derived formula by Ryu et alM. We have seen that we get agreement with the experimental 
results of Majer et ali^ for a value of the anisotropy ~ 525 in good agreement with recent 
measurements of Gaifullin et al. 2 ~ who estimate the anisotropy of their sample to be about 
550. 

There are still many features of BSCCO system in the presence of disorder that are not 
fully understood like the properties of the different phases of the vortex system found in the 
case of point disorder. Future simulations to reproduce the different phases should make 
use of the improved formula for the Josephson interaction in order to get realistic results 
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FIG. 9: Comparison of the results of multilevel Monte Carlo simulations using the new formula for 
the Josephson interaction as compared to the same simulations using the old formula derived by 
Ryu et al. for a magnetic field equal 100 Gauss. The electromagnetic interaction is also included. 
The structure factor (a) and the energy (b) are displayed. 

which are in good agreement with experimental measurements. 
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